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Abstract.  We  consider  propagation  of  waves  in  a  stratified  non-isothermal  magnetic  atmo¬ 
sphere.  The  situation  of  interest  corresponds  to  waves  in  the  outer  solar  (chromosphere  and 
corona)  and  other  stellar  atmospheres.  The  waves  are  simulated  by  using  a  high-resolution, 
well-balanced  finite  volume  based  massively  parallel  code  termed  SURYA.  Numerical  experi¬ 
ments  in  both  two  and  three  space  dimensions  involving  realistic  temperature  distributions, 
driving  forces  and  magnetic  field  configurations  are  described.  Diverse  phenomena  like  mode 
conversion,  wave  acceleration  at  the  transition  layer  and  driving  dependent  wave  dynamics  are 
observed.  We  obtain  evidence  for  the  presence  of  coronal  Alfven  waves  in  some  three  dimen¬ 
sional  configurations.  Although  some  of  the  incident  wave  energy  is  transmitted  into  the  corona, 
a  large  proportion  of  it  is  accumulated  in  the  chromosphere  providing  a  possible  mechanism  for 
chromospheric  heating. 


1.  Introduction 

Waves  and  oscillations  are  a  significant  means  for  the  transport  and  circulation  of  energy  in 
gravitationally  stratified  highly  conducting  astrophysical  plasmas.  Examples  include  waves  emit¬ 
ted  by  localized  sources  within  magnetic  flux  concentrations  such  as  acoustic  sources  in  the  Sun’s 
magnetic  network  and  within  isolated  magnetic  flux  tubes,  knots  and  sunspots.  Other  exam¬ 
ples  pertain  to  waves  in  late  type  stars  and  planetary  magneto- atmospheres.  The  study  of  wave 
propagation  improves  our  understanding  of  the  dynamical  processes  in  the  solar  and  other  stel¬ 
lar  atmospheres  and  contributes  to  explanations  for  phenomena  like  coronal  and  chromospheric 
heating  and  internetwork  oscillations. 

Extensive  studies  of  the  waves  in  the  solar  atmosphere  have  been  performed  with  observations 
from  telescopes  like  the  Swedish  solar  telescope  and  with  satellites  like  SOHO  and  HINODE 
(Banerjee  et  al.  2007).  These  observational  tools  have  also  provided  a  good  deal  of  information 
about  the  detailed  magnetic  structure  of  the  sun.  In  particular,  Van  der  Voort  et  al.  (2005)  show 
that  the  granular  flow  emanating  from  the  solar  convection  zone  arranges  the  magnetic  flux  into 
sheets  that  are  visible  as  thin  bright  features.  These  flux  sheets  are  characterized  by  weak  upflows 
of  the  plasma  inside  them  as  well  as  strong  downflows  in  the  surrounding  medium.  The  flux  sheets 
are  also  subject  to  instabilities.  It  is  well  known  that  convection  results  in  the  excitation  of  acoustic 
waves  at  the  base  of  the  photosphere.  Furthermore,  the  magnetic  instabilities  in  the  photosphere 
excite  magnetic  modes  resulting  in  a  complex  wave  pattern  at  the  base  of  the  photosphere. 

Once  generated  at  the  photospheric  level,  the  waves  travel  up  the  chromosphere  and  into  the 
corona  while  interacting  with  the  magnetic  field  in  a  complicated  manner.  Observational  results 
have  been  obtained  in  McIntosh  et  al.  (2002,  2003,  2004),  Finsterle  et  al.  (2004)  and  references 
therein.  In  McIntosh  et  al.  (2002,  2003)  the  authors  have  found  a  strong  correlation  between 
observations  of  wave  power  from  SOHO  (SUMER)  and  the  magnetic  field  topology  from  SOHO 
(MDI).  The  behavior  of  the  observed  wave  modes  in  the  chromosphere  and  the  corona  has  been 
used  as  a  tool  to  predict  the  magnetic  field  topology.  Such  predictions  were  presented  in  Finsterle 
et  al.  (2004)  for  the  chromosphere  and  in  McIntosh  et  al.  (2004)  for  the  the  corona. 
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Another  interesting  phenomenon  that  the  observations  have  revealed  is  the  existence  of  periodic 
and  spatially  coherent  motions  in  the  photosphere  and  lower  chromosphere.  These  motions  have 
a  peak  in  power  at  around  5  minutes  in  the  photosphere.  Different  5  minute  periodic  motions  in 
the  chromosphere,  transition  region  and  corona  have  been  reported  in  a  series  of  papers,  see  De 
Pontieu  et  al.  (2003a, b),  Finsterle  et  al.  (2008)  and  references  therein. 

Despite  the  wealth  of  observational  data,  it  is  not  yet  possible  to  build  a  consistent  dynamic 
view  of  the  chromosphere  due  to  the  complex  wave  behavior  and  wave- magnetic  field  interaction. 
We  need  to  model  the  conditions  in  the  chromosphere  and  corona  and  use  extensive  numerical 
simulations  to  better  understand  the  observations. 

Theoretical  studies  of  wave  propagation  in  stratified  magneto-atmospheres  have  been  performed 
over  the  last  two  decades.  The  most  popular  modeling  framework  is  based  on  the  equations  of 
magnetohydrodynamics  (MHD)  in  a  stratified  medium.  These  equations  are  a  non-linear  system 
of  PDEs  and  analytical  solutions  are  not  possible,  except  in  the  simplest  of  cases.  In  Zhugzhda  & 
Dzhalilov  ( 1984a, b),  the  authors  obtained  analytical  solutions  for  waves  propagating  in  an  isother¬ 
mal  atmosphere  with  an  oblique  magnetic  field.  More  recent  analytical  studies  were  performed  by 
Cally  (2001,  2006)  by  using  ray  tracing.  These  analytical  studies  are  limited  by  the  simplifying 
assumptions  which  do  not  hold  in  a  realistic  magneto-atmosphere.  Hence,  numerical  simulations 
are  the  key  tools  to  study  wave  propagation  in  the  solar  atmosphere. 

The  design  of  robust  numerical  methods  for  stratified  MHD  equations  is  quite  challenging.  The 
non-linearity  implies  that  solutions  may  contain  discontinuities  in  the  form  of  shock  waves.  The 
shock  structure  of  MHD  is  quite  complicated  as  the  equations  are  neither  strictly  hyperbolic  nor 
convex  (Toth  2000).  Furthermore,  MHD  codes  have  to  handle  the  divergence  constraint  in  a  suit¬ 
able  manner  in  order  to  avoid  numerical  instabilities.  Waves  are  modeled  as  small  perturbations 
of  steady  states  of  the  stratified  MHD  equations.  The  numerical  method  must  be  able  to  cap¬ 
ture  these  perturbations.  Well  balanced  schemes  i.e,  schemes  which  preserve  a  discrete  version  of 
the  steady  state  are  well  suited  for  such  simulations  (Fuchs  et  al.  2010b).  Furthermore,  suitable 
non-reflecting  numerical  boundary  conditions  need  to  be  implemented  as  the  top  boundary  com¬ 
putational  domain  has  to  be  truncated  artificially  at  the  top.  An  elaborate  description  of  these 
difficulties  is  provided  in  Fuchs  et  al.  (2010b,  2011b). 

Numerical  simulations  of  waves  in  the  solar  atmosphere  have  been  presented  in  a  number  of 
papers  in  recent  years.  A  very  partial  list  includes  Rosenthal  et  al.  (2002),  Bogdan  et  al.  (2003), 
Hasan  &  Ballegooijen  (2007,  2008),  Hansteen  et  al  (2006,  2007),  Erdelyi  et  al.  (2007),  Khomenko 
et  al.  (2008),  Fedun  et  al.  (2009),  Carlsson  &  Bogdan  (2006)  and  references  therein.  In  Rosenthal 
et  al.  (2002)  and  Bogdan  et  al.  (2003),  the  authors  study  a  two-dimensional  isothermal  atmosphere 
numerically  using  a  modified  STAGGER  finite  difference  code  of  Nordlund  &  Galsgaard  (1995). 
They  obtain  numerical  evidence  for  the  conversion  of  slow  mode  waves  into  a  combination  of  slow 
and  fast  waves  at  the  magnetic  canopy  i.e  at  plasma  /3  =  1.2.  In  Bogdan  et  al.  (2003),  the  authors 
consider  a  three  dimensional  stratified  isothermal  atmosphere  and  simulate  planar  waves,  using  a 
magnetic  field  extrapolated  from  SOHO/MDI  observations.  The  authors  demonstrate  the  crucial 
role  played  by  the  angle  between  the  magnetic  field  and  the  incident  wave  front  at  the  canopy. 
More  recent  developments  in  this  direction  are  contained  in  Hansteen  et  al.  (2006). 

In  Fedun  et  al.  (2009),  the  authors  simulate  a  three  dimensional  magneto-atmosphere  using  a 
modified  version  of  the  VAC  code  of  Toth  (2000).  The  configuration  includes  both  the  chromo¬ 
sphere  and  the  corona  with  a  temperature  profile  similar  to  the  VAL  IIIC  model  (Vernazza  fetal/ 
1981).  The  background  magnetic  field  is  taken  to  be  a  very  simple  constant  upward  pointing  mag¬ 
netic  field.  One  of  the  key  issues  studied  in  Fedun  et  al.  (2009)  is  the  role  that  driving  frequency 
plays  in  the  wave  propagation. 

Given  that  most  of  the  numerical  studies  have  focused  on  either  simulating  the  chromosphere 
with  a  realistic  magnetic  field  or  on  simulating  the  chromosphere,  transition  region  and  corona  but 
with  a  simplistic  magnetic  field,  we  present  numerical  simulations  for  a  three  dimensional  realistic 
solar  atmosphere  from  the  photosphere  to  the  corona.  The  aim  of  this  paper  is  to  describe  a  code 
for  simulating  the  solar  atmosphere.  This  code  termed  SURYA  has  the  following  properties: 
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(i.)  It  is  a  robust  high  order  finite  volume  scheme  based  on  approximate  Riemann  solvers  of 
the  HLL  type. 

(ii.)  The  code  can  handle  both  two  and  three  dimensional  configurations. 

(iii.)  The  divergence  constraint  is  implemented  through  an  upwind  discretization  of  the  Godunov- 
Powell  source  term  for  ideal  MHD  equations  (Fuchs  et  al.  2011a). 

(iv.)  The  finite  volume  schemes  are  well-balanced  i.e  they  preserve  discrete  equilibrium  states 
for  non-isothermal  atmospheres. 

(v.)  Realistic  temperature  profiles,  background  magnetic  fields  and  photospheric  driving  are 
incorporated. 

(vi.)  Neumann  type  numerical  boundary  conditions  (Fuchs  et  al.  2010a)  are  employed  at  the 
top  boundary. 

(vii.)  The  code  is  implemented  on  massively  parallel  architectures. 

Many  features  of  the  code  are  novel.  In  particular,  the  code  is  based  on  the  first  well-balanced 
schemes  for  stratified  MHD  equations  (Fuchs  et  al.  2010b,  2011b). 

We  perform  and  report  several  numerical  simulations  with  SURYA  in  this  paper.  The  simula¬ 
tions  are  based  on 

•  A  non-isothermal  stratified  MHD  model  that  includes  a  synthetic  (inspired  by  the  VAL- 
IIIC)  temperature  profile  for  the  chromosphere  and  the  corona. 

•  Synthetic  magnetic  fields  (following  Bogdan  et  al.  2003)  in  both  two  and  three  space 
dimensions. 

•  A  magnetic  field  extrapolated  from  SOHO/MDI  that  can  replace  the  three-dimensional 
synthetic  field. 

•  Photospheric  driving  mechanisms  of  different  types  and  with  different  frequencies. 

We  focus  on  how  the  waves  induced  at  the  photosphere  travel  up  the  chromosphere  and  the 
corona  and  how  they  interact  with  the  underlying  magnetic  field.  In  particular,  we  focus  on  how 
much  energy  these  waves  carry  up  the  atmosphere  and  deposit  in  the  chromospheric  and  coronal 
plasmas.  A  related  issue  that  we  address  is  the  role  played  by  the  frequency  of  the  driving  in  the 
energy  transfer.  Most  of  the  analysis  is  performed  in  two  space  dimensions  and  with  synthetic 
magnetic  fields  in  three  dimensions.  We  also  include  a  simulation  with  an  extrapolated  magnetic 
field  from  SOHO/MDI  observations  and  with  photospheric  driving,  also  from  SOHO  data,  in  order 
to  demonstrate  the  ability  of  our  code  to  handle  realistic  scenarios. 

The  rest  of  this  paper  is  organized  as  follows:  in  section  [2j  we  present  the  model  equations  and 
the  steady  states  of  interest.  The  numerical  schemes  are  described  in  section [3] and  in  Appendix  A. 
Numerical  simulations  in  two  and  three  dimensions  are  presented  in  sections  [4]  and  [5|  respectively. 


2.  The  Model 

We  model  the  outer  solar  atmosphere  with  a  modified  version  of  the  stratified  MHD  equations 
(Fuchs  et  al.  2010b).  In  non-dimensionalized  form,  the  governing  equations  are 

pt  +  div(pu)  =  0, 

(pn)t  +  div  (pu  0  u  +  (p  +  t  |B|2  +  B  •  B^)  I  —  B(g>B  —  B  0  B  —  B  0  B^) 

=  -  (b  +  b)  (div  B)  -  pge3, 

(2-1)  .  „  „  ,  V  7 

B t  +  div  fu(g)B  —  B(g)u  +  u(g)B  —  B(g)uj  =  — u(divB), 

Et  +  div  (Je  +  p  +  t  |B|2  +  B  •  u  -  (u  •  B)B  -  (u  •  fi)  B^ 

=  —  (u  •  B)(div  B)  -  pg( u  •  e3) . 

Here,  p  is  the  density,  u  =  (ui,U2,us)  the  velocity,  p  the  thermal  pressure,  g  the  constant  accel¬ 
eration  due  to  gravity,  and  e3  represents  the  unit  vector  in  the  vertical  (z-)  direction.  The  total 
magnetic  field  is  B  +  B,  where  B  =  (B i,  F?2,  H3)  is  the  deviation  from  a  background  potential  field 
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B,  i.e.  B  satisfies  the  following  assumptions, 

(2.2)  B t  =  0,  div(B)  =  0,  and  curl(B)  =  0. 

The  total  energy  E  is  determined  by  the  ideal  gas  equation  of  state: 


E  = 


V 


W  +  iplul2. 


7-1  2 

Writing  <ED  as  a  balance  law  explicitly,  we  obtain, 

(2.3)  Ut  +  (f(U,  B)*  +  g(U,  B)„  +  h(U,  B)z  =  s^U,  B)  +  s2(U,  B)  +  +s3(U,  B)  +  s9(U), 
where 

U  =  {p,  pui ,  pu2 ,  pus  ,B1,B2,Bs,E} 

is  the  vector  of  conserved  variables  and  B  =  {£?i,  B2,  B^}  is  the  background  magnetic  field.  The 
fluxes  f ,  g,  h,  Godunov- Powell  source  terms  s1,2,3  and  the  gravity  source  term  can  be  read  from 

& 


2.1.  Steady  states.  We  construct  equilibrium  solutions  of  (2.1 )  by  setting  the  velocity  u  and  the 
magnetic  field  deviation  B  to  zero.  From  the  ideal  gas  equation  of  state,  we  obtain 

(2.4)  p  =  gHpT , 

for  a  constant  H  and  the  temperature  T.  We  assume  that  the  steady  state  temperature  T  =  T(z) 
varies  only  in  the  vertical  direction.  Substituting  (2.4)  in  (2.1)  and  assuming  u  =  0  leads  to 


(2.5) 


dp 

dz 


V 


HT(z)' 


_  ol(z) 


Solving  the  above  equation  explicitly  yields 

p(x,  y,  Z )  =  p(z)  =  p0e  . 

Here,  po  is  a  constant  and 

rz  1 

(2.6)  a(x,y,z)  =  a(z)=  Tf^ds. 

Jo  T{s) 

Similarly,  we  can  calculate  the  steady  state  density  as 

p{x,y,z)  =  P\Zj  =  ——e  m  , 

1  [z) 

with  (po,To)  being  constant.  Combining  the  above  expressions,  we  obtain  the  steady  state: 

(2.7) 


_  /  \  O0T0  a(z)  a(z) 

u  =  0,  B  =  0  p(z)  =  e  H  ,  p(z)  =  poe  H  . 


T(z ) 

The  steady  state  pressure  and  density  are  scaled  in  terms  of  the  function  a  which  in  turn  depends 
on  the  temperature.  Furthermore  cu  is  a  monotonically  increasing  function  as  the  temperature  is 
always  positive.  In  the  simplest  case  of  an  isothermal  atmosphere,  i.e.,  T  =  T  for  some  constant 
T,  the  expression  (2.7)  simplifies  as 

(2.8)  u  =  0,  B  =  0,  p(x,y,z)  =  p0e~™,  p(x,y,  z)  =  p0e~™  , 

and  the  pressure  and  density  decay  exponentially. 


2.1.1.  Hydrodynamic  steady  state.  We  obtain  a  purely  hydrodynamic  steady  state  by  assuming 
that  the  background  magnetic  field  is 

(2.9)  B  =  0. 
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2.1.2.  Magnetic  steady  states.  Non-trivial  solutions  of  (2.2)  lead  to  interesting  magnetic  steady 
states.  Note  that  solutions  of  (2.2)  can  be  characterized  by  vector  harmonic  functions. 

In  three  space  dimensions,  a  Fourier  solution  of  (2.2)  can  be  written  down  explicitly  as, 

(2.10) 

L  M  ] 

'  alm  ^-2t TVl2+m2z 


Bi(x,y,z)  = 


1=1  m= 0 


L  M 


B2(x,y,z)  = 


Vl2  +  rn2 


m  aim  ^-2tt yi2+m2z  I 


7  n  ,  yjl2  +  - 

1=0  m= 1 


L  M 


Bz(x,  y,  z)  = 


EE  CUmT 

1=0  m= 0 


—  27 x\/l2-\-m2  z 


aim  sin(27r lx)  cos(27r my)  —  bim  cos(27 xlx)  sin(27r my) 
cim  cos(27 xlx)  cos(27rrm/)  +  dim  sin(27r/x)  sin(27r my)^j , 
aim  cos(27 xlx)  sin(27rm7/)  —  bim  sin(27r lx)  cos(2imiy) 

+  cim  sin(27r/x)  sin(27rrra/)  —  dim  cos(27 xlx)  cos(27rrm/)^ , 
^  aim  cos(27 xlx)  cos(27t my)  +  bim  sin(27r lx)  sin(27r my) 

+  cim  sin(27r lx)  cos(27t my)  +  dim  cos(27 xlx)  sin(27r m?/)J . 


Here,  a/m,  Qm,  are  the  Fourier  coefficients  corresponding  to  the  background  magnetic  field 
Bs(x,y,0)  at  the  bottom  of  the  domain,  and  L,  M  are  the  maximum  number  of  modes  for  the 
indices  l  and  m  respectively.  The  factor  cpm  is  1/4  if  Z  =  m  =  0,  1/2  if  l  or  m  is  zero,  and  1 
otherwise.  It  can  be  readily  checked  that  (2.10)  satisfies  (2.2).  See  also  Fuchs  et  al.  (2010a, b)  for 
a  two  dimensional  version. 


3.  The  code 


We  approximate  (2.1)  with  a  second-order  accurate  finite  volume  scheme.  The  scheme  has 


the  added  advantage  of  being  well  balanced  i.e,  it  preserves  a  discrete  version  of  the  steady  state 
(2.7)  for  any  background  magnetic  field  B  satisfying  (2.2)  and  for  any  steady  state  temperature 


distribution  T  =  T(z).  The  scheme  is  described  in  detail  in  Appendix  A. 

This  well-balanced  high-order  scheme  is  implemented  in  a  modular  C++  based  code  termed 
The  code  includes  a  set  of  approximate  Riemann  solvers,  high-order  non-oscillatory 
reconstruction  and  time  integration  routines.  Realistic  initial  and  boundary  conditions  are  also 
specified.  A  wide  range  of  background  magnetic  fields  B  are  included.  The  code  is  parallelized  with 
the  MPI  library,  using  a  domain  decomposition  technique.  The  parallelization  is  straightforward 
as  the  schemes  are  explicit  and  do  not  need  any  staggered  grids.  A  python  front  end  to  the  code 
is  included  for  configuring  data  and  results  are  visualized  using  matplotlib  for  two-dimensional 
simulations  and  MAYAVI2  for  three-dimensional  simulations.  All  the  results  presented  below  are 
from  experiments  performed  on  the  TITAN  cluster  of  the  University  of  Oslo  and  on  the  STALLO 
cluster  of  the  University  of  Tromsp. 


4.  Numerical  results  in  two  space  dimensions 

All  the  quantities  that  appear  in  are  non-dimensionalized  suitably  from  the  realistic  solar 
parameters  used  in  Bogdan  et  al.  (2003),  shown  in  table  [lj  The  constants  are  acceleration  due 
to  gravity,  g  =  2.74,  the  thermodynamical  parameter  H  =  0.158  and  bottom  initial  pressure 
Po  =  1.13.  All  subsequent  two-dimensional  experiments  are  performed  on  the  domain  [x,z]  G 
[0,4]  x  [0,8], 

4.1.  Temperature  profile:  We  use  a  temperature  profile  for  the  chromosphere  and  corona  that  is 
based  on  the  VAL-IIIC  model  (Vernazza  et  al.  1981),  see  figure [l] (left  panel).  We  fit  the  following 


1  http :  /  /  folk,  uio .  no  /  mcmurry  /  amhd  / 
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Table  1.  Units  used  in  this  paper 


Quantity  CGS  units 


length 
time 
velocity 
pressure 
density 
magnetic  field 


108  cm 

102  s 
106  cm  s_1 
105  g  cm-1  s-2 
10-11  g  cm-3 
1120  G 


(simplified)  temperature  profile, 


1,  if  z  <  1 

I  +792 (z  -  l)2  +  1,  if  1  <  *  <  1.25 

^  ~  |  -792 (z  -  1.5)2  +  100,  if  1.25  <  z  <  1.5 

,100,  if  1.5  <  z 

to  the  VAL-IIIC  profile.  The  domain  starts  from  the  photosphere  and  continues  through  the 
chromosphere  into  the  corona.  The  transition  region  is  between  z  =  1  and  z  =  1.5.  The  equations 
solved  in  our  code  will  not  accurately  represent  coronal  physics.  In  this  paper  we  are  interested  in 
the  waves  that  are  transmitted  through  the  transition  region,  rather  than  in  coronal  dynamics. 


temperature[5000  K]  density  [g/cm^B] 


(a)  Real  sun  (VAL  IIIC  profile) 


(b)  Model  sun 


Figure  1.  Steady  state  temperature  distribution  for  the  real  and  model  solar  atmospheres. 


4.2.  Background  magnetic  field.  In  two  space  dimensions,  we  choose  a  realistic  two-dimensional 


background  magnetic  field  of  the  form  (2.10)  as  follows:  We  let  B3(x,  0,0)  approximate 


(4.1) 


B3(x,  0,0)  =  2.7e 


—  7.2+ 


1.3e-4°(r-0'6)2,  r  =  \x-2\,x  e  [0,4] 


up  to  the  first  fourteen  terms  in  the  Fourier  series.  The  full  magnetic  field  B (x,y,  z)  follows  from 
the  potential  field  assumption  as  in  (2.10).  The  resulting  potential  field  consists  of  a  large  unipolar 


magnetic  flux  concentration  surrounded  on  each  side  by  two  smaller  concentrations  of  opposite 
polarity  field  (Bogdan  et  al.  2003,  figure  1  for  an  illustration).  The  initial  conditions  are  set  to 
the  steady  state  (|2.7|).  We  consider  the  following  configurations: 


4.3.  High  frequency  planar  waves  with  vertical  forcing.  With  the  above  background  tem¬ 
perature  profile  and  magnetic  field,  we  excite  planar  waves  at  the  bottom  boundary  by  prescribing 

(4.2)  u(x,0 ,  £)  =  (0, 0,  0.1  sin(27rco’t))T. 
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Here,  the  frequency  of  the  driving  is  the  parameter  uj.  For  our  first  experiment,  we  choose  uo  =  3.0. 
In  figure  [2j  we  plot  time  snapshots  (at  different  time  levels)  of  the  velocity  component  parallel  to 
the  magnetic  field  and  the  velocity  component  perpendicular  to  the  magnetic  field.  The  plots  show 
the  planar  waves,  excited  at  the  photosphere,  traveling  up  the  chromosphere  and  into  the  corona. 
The  velocity  in  the  direction  of  the  magnetic  field  is  an  indicator  of  both  slow  and  fast  mode  waves 
(Bogdan  et  al.  2003)  whereas  the  velocity  perpendicular  to  the  magnetic  field  detects  only  fast 
mode  waves.  As  shown  in  the  figure,  we  excite  acoustic  type  planar  waves  at  the  bottom  boundary. 
This  is  expected  as  the  vertically  directed  forcing  has  the  same  direction  as  the  magnetic  field  for 
most  of  the  bottom  boundary.  There  are  small  regions  near  the  center  (symmetric  about  x  =  2) 
where  the  magnetic  field  turns  and  is  not  aligned  to  the  vertical  direction,  causing  some  disruption 
of  the  largely  planar  wave  front  structure.  Clearly,  the  velocity  in  the  direction  of  the  magnetic 
field  is  symmetric  whereas  the  one  in  the  transverse  direction  is  antisymmetric.  Furthermore,  the 
velocity  in  the  transverse  direction  is  strong  exactly  where  the  velocity  in  the  direction  of  the 
magnetic  field  is  weak. 

The  waves  in  the  direction  of  the  magnetic  field  are  pushed  up  vertically  whereas  the  ones  in 
the  transverse  direction  are  pushed  up  sideways.  Once  the  initial  wave  hits  the  magnetic  canopy 
(/ 3  =  1.2),  the  initial  acoustic  wave  mode  is  converted  into  a  combination  of  fast  and  slow  mode 
waves.  Note  that  this  mode  conversion  is  taking  place  in  the  (horizontal)  center  of  the  domain 
as  the  magnetic  field  is  strongest  here.  The  resulting  waves  have  an  elongated  structure  as  the 
Alfven  speed  is  greater  than  the  sound  speed  in  this  region. 

Then,  the  waves  hit  the  transition  layer  and  are  accelerated  due  to  the  increase  in  sound  speed 
by  an  order  of  magnitude  as  the  temperature  rises  by  two  orders  of  magnitude.  In  particular,  the 
fast  waves  zoom  out  of  the  transition  layer  with  the  wave  signature  visible  in  the  perpendicular 
component  of  velocity  at  time  £  =  0.8.  At  £  =  1.2,  the  velocity  in  the  direction  of  the  magnetic 
field  has  a  complex  structure  with  rich  spatial  variation.  The  magnetic  field  along  the  sides,  i.e. 
in  the  region  of  [0, 1]  x  [2,  8]  U  [3, 4]  x  [2,  8]  is  almost  vertical  throughout  the  model,  and  has  little 
effect  on  the  shape  of  the  waves.  The  region  where  the  acoustic  type  waves  predominate  and  the 
region  with  a  more  complex  wave  structure  are  clearly  delimitated  in  figure  |2jc)  by  the  j3  m  0.1 
isoline.  The  more  complex  magnetic  field  at  the  base  in  x  G  [1,3]  results  in  the  more  complex  wave 
structure  in  this  central  region.  The  fast  mode  waves  are  visible  in  the  perpendicular  component 
of  velocity  and  have  a  compressive  saw-tooth  like  behavior.  This  compression  is  also  visible  in  the 
normal  component  indicating  the  formation  of  shock  waves.  Furthermore,  the  waves  generated  at 
the  center  are  spreading  out.  Being  fast  waves,  they  can  travel  across  magnetic  field  lines. 

As  we  are  interested  in  how  the  waves  transport  energy  from  the  photosphere  into  the  upper 
chromosphere  and  the  corona,  we  consider  the  relative  change  in  total  energy  (kinetic  +  internal 
+  magnetic  energy),  i.e 

cR  (4-\ _ Ejj(t)  —  Ejj( 0) 

iA  }  ■  Eitj( 0)  • 

The  relative  change  is  a  good  indicator  of  how  energy  is  transported.  Since  we  are  interested 
in  studying  the  spatial  variation  in  the  energy  transfer,  the  total  energy  at  different  horizontal 
transects  is  plotted  vs.  time  in  figure  [3J  We  choose  four  different  horizontal  locations,  namely 
x  =  0.7,  2.0,  2.6,  3.8.  Figure  [3]  shows  a  rich  spatial  variation  in  the  change  of  energy.  At  the 
locations  x  =  0.7  and  x  =  3.8,  the  basal  magnetic  field  is  vertical  and  quite  weak.  The  energy 
transport  is  similar  in  the  two  locations.  In  particular,  the  original  acoustic  wave  travels  up  the 
chromosphere  and  hits  the  transition  region  without  much  interaction  with  the  magnetic  canopy. 
At  the  bottom  of  the  transition  layer,  the  wave  undergoes  both  refraction  and  reflection.  A  part 
of  the  energy  is  transmitted  up  into  the  corona.  The  refracted  wave  has  a  higher  speed  due  to  the 
larger  temperature  in  the  corona.  A  part  of  the  wave  energy  is  also  reflected  from  the  transition 
region  and  travels  down.  This  wave  interacts  with  the  second  upward  moving  incident  wave  and 
forms  a  complex  interference  pattern.  Another  part  of  the  wave  energy  from  the  first  incident  wave 
moves  horizontally  along  the  transition  layer.  The  second  and  third  incoming  waves  have  similar 
behavior  but  are  affected  by  the  waves  reflected  from  the  transition  layer  resulting  in  interference 
patterns. 
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<U,B>/|B|  t=0.4 


<U,B>/|B|  t=0.8 


(b)  T  =  0.8 


<U,B>/|B|  t=1.2 


(c)  T=  1.2 


(d)  T  =  0.4 


. 


0.0 


(f)  T=  1.2 


Figure  2.  Results  for  the  high  frequency  planar  wave  experiment.  Top:  Velocity 
in  the  direction  of  the  magnetic  field.  Bottom:  Velocity  perpendicular  to  the 
direction  of  the  magnetic  field.  All  the  results  are  with  a  second-order  accurate 
scheme  on  a  mesh  of  400  x  800  points. 


At  the  location  x  =  2,  the  magnetic  field  is  strong  at  the  bottom.  Here,  we  see  a  pronounced 
difference  in  wave  behavior  from  the  previous  horizontal  locations.  In  particular,  the  first  incident 
slow  mode  wave  hits  the  magnetic  canopy  and  generates  fast  waves.  The  fast  mode  wave  is  clearly 
visible  at  t  =  0.7.  It  reaches  the  transition  layer  and  is  accelerated.  This  wave  carries  only  a  very 
small  fraction  of  the  wave  energy.  The  slow  mode  wave  hits  the  transition  region  and  a  part  of  it 
is  transmitted  (refracted)  into  the  corona.  A  part  of  the  wave  energy  is  also  reflected  downwards 
and  some  of  the  incident  wave  energy  remains  within  the  transition  layer.  Notice  that  the  reflected 
wave  interacts  with  the  fast  mode  excited  from  the  second  incident  wave  at  time  t  —  1.1.  The  waves 
excited  later  from  the  bottom  boundary  show  similar  behavior.  In  this  case,  the  reflection  from 
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(a)  x  =  0.7 


(b)  x  =  2.0 


(c)  x  =  2.6 


(d)  x  =  3.8 


Figure  3.  Results  for  the  high  frequency  planar  wave  experiment.  The  relative 
change  in  total  energy  as  a  function  of  time  and  height  is  plotted  at  different 
horizontal  locations.  All  the  results  are  with  a  second-order  accurate  scheme  on 
a  mesh  of  400  x  800  points. 


the  transition  layer  appears  to  be  of  smaller  amplitude  than  the  previous  horizontal  locations.  At 
the  horizontal  location  x  =  2.6,  the  magnetic  field  is  of  moderate  strength  at  the  bottom  and  the 
waves  transport  energy  in  a  similar  manner  to  the  x  =  2  case.  Since  the  waves  hit  the  canopy  at  an 
angle  to  the  magnetic  field,  we  should  expect  significant  mode  conversion,  Indeed,  fast  mode  waves 
are  clearly  excited  at  the  magnetic  canopy.  Waves  are  transmitted  through  the  transition  layer 
into  the  corona.  However,  a  considerable  amount  of  the  wave  energy  remains  within  the  transition 
region.  We  observe  particularly  complicated  interference  patterns  here,  which  can  be  attributed 
to  the  presence  of  both  mode  conversion  at  the  canopy  and  reflections  from  the  transition  region. 

Although  the  waves  show  rich  spatial  variation,  there  are  a  some  common  features.  In  particular, 
the  waves  carry  some  of  the  incident  energy  into  the  corona.  However,  most  of  the  energy  is  dumped 
at  the  base  of  the  transition  region  and  in  the  upper  chromosphere.  This  is  clearly  seen  in  figure 
[4]  where  we  plot  the  relative  change  in  energy,  integrated  across  all  horizontal  locations,  in  time. 
This  deposition  of  energy  at  the  base  of  the  transition  layer  and  the  chromosphere  is  a  plausible 
candidate  for  explaining  chromospheric  heating. 
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E,  f  =  3.0,  angle  at  tp  =  47 


0.2  0.4  0.6  0.8  1.0  1.2  1.4  1.6 

Figure  4.  Relative  change  in  total  energy,  integrated  across  all  horizontal  loca¬ 
tions,  vs.  time  in  the  high  frequency  vertical  driving  case. 


4.4.  Low  frequency  planar  waves.  We  change  the  driving  frequency  in  (4.2)  to  uo  =  0.3.  Thus, 
we  drive  waves  at  the  photosphere  with  a  frequency  which  is  an  order  of  magnitude  less  than  the 
previous  experiment.  The  velocities,  in  directions  parallel  and  perpendicular  to  the  magnetic  field 
are  plotted  in  figure  [5]  The  results  are  quite  similar  to  the  high  frequency  case  (see  figure  [5]  for  a 
comparison).  Note  that  the  time  levels  for  the  snapshots  in  the  two  figures  are  different. 

A  key  difference  with  the  high  frequency  excitation  lies  in  how  the  energy  is  trapped  by  the 
transition  region.  In  order  to  illustrate  this  process,  we  plot  the  relative  change  in  energy  for  this 
case  and  compare  with  the  high  frequency  excitation  in  figure  [6j  As  shown  in  this  figure,  there 
are  major  differences  between  both  cases.  In  particular,  there  is  a  greater  buildup  of  the  energy  at 
the  base  of  the  transition  region  and  in  the  chromosphere  for  low  frequency  waves  than  for  high 
frequency  waves.  It  is  not  unexpected  that  more  wave  energy  is  trapped,  since  the  frequency  of  the 
forcing  at  the  bottom  is  below  the  acoustic  cut-off  of  the  photosphere.  Furthermore,  this  buildup 
is  spread  over  a  larger  portion  of  the  base  of  the  transition  region  in  the  low  frequency  case.  This 
suggests  that  there  will  be  greater  chromospheric  heating  with  low  frequency  excitations  than  with 
high  frequency  ones. 


4.5.  Transverse  forcing.  In  the  above  numerical  experiments,  the  forcing  was  in  the  normal 
direction.  In  this  experiment,  we  prescribe  the  transverse  driving  force, 

(4.3)  u(x,0 ,£)  =  (0.1  sin(27ra;t),  0,  0)T. 

Here,  u  is  the  driving  frequency.  We  start  by  considering  high  frequency  excitations  at  uj  =  3.0. 
The  numerical  results  for  this  experiment  are  presented  in  figures  |7|8|  and  [9]  In  figure  [7J  we 
plot  the  velocities  in  the  direction  of  the  magnetic  field  and  perpendicular  to  the  magnetic  field. 
There  is  a  clear  difference  in  the  wave  behavior  due  to  transverse  driving  from  the  vertical  driving 
case  above.  Figure  [7]  should  be  compared  with  figure  [2]  We  excite  a  combination  of  fast  and 
slow  mode  waves  at  the  bottom  boundary.  The  waves  in  the  direction  of  the  magnetic  field  are 
anti-symmetric  about  the  x  =  2  line,  whereas  those  in  the  direction  perpendicular  to  the  field  are 
symmetric  about  the  x  =  2  line.  Unlike  the  vertical  driving,  we  see  that  the  waves  behave  as  if 
they  are  generated  by  a  localized  piston,  located  at  the  middle  of  the  bottom  boundary.  There  is 
very  little  wave  excitation  in  the  region  where  the  magnetic  field  is  weak  at  the  bottom,  although 
the  forcing  acts  all  across  the  bottom  boundary.  The  waves  move  up  the  chromosphere  and  hit 
the  magnetic  canopy  where  mode  conversion  takes  place.  The  waves  are  accelerated  when  they  hit 
the  transition  layer.  The  resulting  waves  have  a  semi-circular  arch  shape  and  have  a  pronounced 
sideways  turning  motion,  which  should  be  contrasted  with  the  vertical  driving  case. 
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<(U1,U2),(-B2,B1)>/|B|  t=0.8 


(d)  T  =  0.8 


<(U1,U2),(-B2,B1)>/|B|  t=2.4 


(f)  T  =  2.4 


Figure  5.  Results  for  the  low  frequency  planar  wave  experiment.  Top:  Velocity 
in  the  direction  of  the  magnetic  field.  Bottom:  Velocity  perpendicular  to  the 
direction  of  the  magnetic  field.  All  the  results  are  with  a  second-order  accurate 
scheme  on  a  mesh  of  400  x  800  points. 


The  relative  energy  change  (as  a  function  of  time  and  height)  at  four  different  horizontal  loca¬ 
tions  x  =  0.7,  2.0.2,  6,  3.8  is  shown  in  figure  [8]  Again,  there  are  considerable  differences  between 
this  example  and  the  one  with  vertical  driving.  In  particular,  the  amplitude  of  energy  gains  at  all 
four  locations  is  considerably  smaller  in  this  case  (compare  with  figure  [3]).  Furthermore,  a  much 
larger  proportion  of  the  energy  leaks  into  the  corona.  There  is  again  a  rich  spatial  variation  in 
the  energy  transfer.  At  x  =  2,  the  magnetic  field  is  strong  at  the  bottom  and  the  initial  waves  are 
acoustic  type  waves.  Fast  mode  waves  are  also  generated  when  the  incident  waves  hit  the  magnetic 
canopy.  These  fast  waves  are  transmitted  into  the  corona.  As  the  slow  wave  hits  the  transition 
region,  a  part  of  the  energy  is  reflected.  Notice  that  the  reflections  have  lower  amplitude  in  this 
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relchange  E  relchange  E 


(a)  u>  =  3.0,  T  =  1.2 


(b)  u  =  0.3,  T  =  2.4 


Figure  6.  A  comparison  of  the  relative  change  in  energy  for  the  high  frequency 
case  and  the  low  frequency  case.  Both  results  are  for  the  vertical  driving  at 
comparable  times. 


case  than  the  one  with  vertical  forcing.  Furthermore,  some  of  the  wave  energy  remains  within 
the  transition  layer.  The  trailing  waves  interact  with  the  reflected  waves  and  horizontal  waves  in 
the  transition  layer,  forming  a  complex  pattern.  A  similar  qualitative  behavior  is  encountered  at 
the  spatial  location  x  =  2.6.  The  main  difference  with  the  previous  location  lies  in  the  fact  that 
fewer  fast  waves  are  excited.  At  x  =  0.7,  the  magnetic  field  is  approximately  zero  at  the  bottom 
z  =  0.  There  are  few  visible  fast  waves  and  most  of  the  energy  is  carried  by  slow  mode  waves.  A 
large  fraction  of  the  waves  seen  around  the  /3  =  1.2  isoline  after  time  t  ss  1.2,  are  waves  generated 
around  the  center,  entering  form  the  side,  compare  figure  [lOj  They  interact  with  upward  moving 
waves  forming  a  complex  interference  pattern.  At  the  location  x  =  3.8,  the  magnetic  field  is 
again  close  to  zero  at  the  bottom  boundary  and  we  see  a  similar  behavior  to  the  case  x  =  0.7.  A 
large  proportion  of  the  incident  wave  energy  is  transmitted  into  the  corona.  The  relative  change, 
integrated  in  the  horizontal  direction,  (as  a  function  of  time)  is  shown  in  figure  [9j  It  highlights 
the  differences  with  the  vertical  driving  case  (figure  [4]) .  The  amplitude  of  energy  change  is  smaller 
in  this  case  but  a  larger  proportion  of  the  wave  energy  leaks  into  the  corona. 

We  repeated  the  simulations  with  the  driving  (|4.3|)  but  with  frequency  uj  =  0.3.  The  wave 
dynamics  for  this  low  frequency  case  was  very  similar  to  the  one  with  high  frequency  excitations. 
The  main  difference  in  the  relative  change  in  energy  is  shown  in  figure  [lO]  As  in  the  vertical 
driving  case,  the  low  frequency  excitations  result  in  a  greater  buildup  of  energy  at  the  base  of  the 
transition  region  and  in  the  chromosphere. 
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<U,B>/|B|  t=0.4 


<U,B>/|B|  t=0.8 


<U,B>/|B|  t=1.2 


(a)  T  =  0.4 


(b)  T  =  0.8 


(c)  T  =  1.2 


(d)  T  =  0.4 


(e)  T  =  0.8 


(f)  T  =  1.2 


Figure  7.  Results  for  the  high  frequency  planar  wave  experiment  with  driving 


(4.3).  Top:  Velocity  in  the  direction  of  the  magnetic  field.  Bottom:  Velocity 
perpendicular  to  the  direction  of  the  magnetic  field.  All  the  results  are  with  a 
second-order  accurate  scheme  on  a  mesh  of  400  x  800  points. 


5.  Numerical  results  in  three  space  dimensions 

In  this  section,  we  will  present  numerical  results  for  wave  propagation  in  a  three  dimensional 
model  solar  atmosphere.  The  temperature  profile  is  exactly  the  same  as  in  the  previous  section. 

5.1.  Synthetic  magnetic  field.  To  begin  with,  we  consider  a  synthetic  background  magnetic 
field  of  the  following  form.  Let 

(5-1) 

B3{x,y,  0)  =  0.0275  (2.8e“6-4r 2  -  0.7e-40(r-°-75)2)  ,  x,y  e  [0,4],r  =  y/(x  -  2)2  +  (y  -  2)2. 

In  3d  we  take  the  zeroth  component  of  the  Fourier  term  (2.10|)  4  times,  i.e.  do,o  =  4ao,o-  The 
magnetic  field  is  constructed  from  the  first  16  x  16  Fourier  co-efficients.  The  computational  domain 
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E  at  x  =  2.6,  f  =  3.0,  angle  at  tp  =  35  E  at  x  =  3.8,  f  =  3.0,  angle  at  tp  =  4 


(c)  X  =  2.6 


(d)  x  =  3.8 


Figure  8.  Results  for  the  high  frequency  planar  wave  experiment  with  transverse 
driving  (4.3).  The  relative  change  in  total  energy  as  a  function  of  time  and  height 
is  plotted  at  different  horizontal  locations,  characterized  by  different  attack  angles. 
All  the  results  are  with  a  second-order  accurate  scheme  on  a  mesh  of  400  x  800 
points. 


E,  f  =  3.0,  angle  at  tp  =  16 


Figure  9.  Relative  change  in  total  energy,  integrated  across  all  horizontal  loca¬ 
tions,  vs.  time  in  the  2-D  transverse  driving  case. 
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(a)  lj  =  3.0,  T  =  1.2  (b)  cj  =  0.3,  T  =  2.4 

Figure  10.  A  comparison  of  the  relative  change  in  energy  for  the  the  high  fre¬ 
quency  case  and  the  low  frequency  case.  Both  results  are  for  transverse  driving 
at  comparable  times. 

is  (x,  y,  z )  e  [0, 4]  x  [0, 4]  x  [0,8]  and  the  parameters  g,  H  and  p0  are  the  same  as  in  the  previous 
section.  Figure  [TT]  shows  the  field  lines  of  the  background  magnetic  field  in  white  and  the  magnetic 
field  strength  at  z  =  0.  Furthermore,  the  blue  isosurface  is  the  location  of  the  magnetic  canopy, 
i.e.  (3  =  1.2. 

5.1.1.  Planar  waves  with  vertical  forcing.  We  excite  planar  waves  at  the  bottom  boundary  with 
the  forcing, 

(5.2)  u(x,  y,  0,  t)  =  (0,  0,  0.1  sin(27rcc;t))T, 

with  the  frequency  uj  set  to  3.  The  results  of  this  experiment  are  presented  in  figures  [121  [13]  and 


in  the  direction  of  the  magnetic  field  as  well  as  in  two  perpendicular  transverse  directions  to  the 
magnetic  field.  These  two  directions  are  calculated  as  in  Carlsson  and  Bogdan  (2006)  and  are 
termed  as  the  principal  normal  direction  and  the  binormal  direction.  This  decomposition  of  the 
velocity  field  serves  to  investigate  the  existence  of  two  transverse  wave  modes. 

We  show  isosurfaces  of  the  velocity  field  in  all  the  three  directions  in  figure  [12]  The  figure 
shows  that  there  is  a  marked  difference  between  velocity  in  the  normal  direction  and  in  the  two 
transverse  directions.  In  particular,  the  velocity  in  the  normal  direction  shows  waves  that  are 
excited  all  over  the  domain,  including  at  the  edges  where  the  magnetic  field  is  quite  weak  at  the 
bottom.  On  the  other  hand,  the  transverse  waves  are  present  only  in  the  center  of  the  domain 
where  the  basal  magnetic  field  is  strong.  This  is  similar  to  the  two  dimensional  situation.  The  two 
transverse  wave  modes  are  quite  different  from  each  other.  The  mode  illustrated  by  the  principal 
normal  direction  has  a  four  fold  (quadrant)  symmetry  whereas  the  waves  in  the  binormal  direction 


14  As  there  are  three  directions,  we  follow  Carlsson  and  Bogdan  (2006)  and  show  the  velocity 
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Figure  11.  Background  magnetic  field  in  3  dimensions.  The  magnetic  field  lines 
are  shown  in  white.  At  z  —  0  we  show  the  strength  of  the  magnetic  field.  The 
points  where  /5  =  1.2  he  on  the  blue  isosurface. 


have  a  (approximately)  radial  symmetry.  These  differences  in  the  two  transverse  directions  may 
indicate  the  existence  of  two  different  transverse  wave  modes  (Carlsson  &  Bogdan  2006).  We  will 
examine  this  issue  in  some  detail  later  in  this  section. 

The  relative  change  in  total  energy  as  function  of  time  and  at  four  different  horizontal  locations 
is  shown  in  figure  [13]  The  figures  show  pronounced  spatial  variation  in  how  the  waves  transport 
energy  up  the  solar  atmosphere.  At  the  location  (2.0,  2.0),  i.e  the  center  of  the  horizontal  domain, 
the  magnetic  field  is  strong  at  the  bottom.  Here,  we  see  that  the  initial  acoustic  waves  carry  energy 
up  the  chromosphere.  The  magnetic  field  is  strong  enough  that  the  waves  do  not  intersect  the 
iso-surface  (3  =  1.2  at  this  location.  Consequently,  there  are  no  visible  signs  of  mode  conversion. 
The  waves  hit  the  transition  layer  and  a  part  of  the  wave  energy  is  transmitted  (refracted)  into  the 
corona.  These  waves  are  accelerated  at  the  transition  region.  There  is  little  evidence  of  reflections 
from  the  transition  region  in  the  figure.  This  phenomenon  could  well  be  explained  by  the  reflected 
waves  traveling  along  some  other  direction  and  therefore  not  being  visible  at  the  location  (2.0,  2.0). 
Still,  a  large  proportion  of  the  wave  energy  is  dumped  at  the  base  of  the  transition  region.  Very 
similar  wave  dynamics  is  observed  at  the  locations  (2.25,  2.0)  and  (2.2,  2.2)  where  the  magnetic 
field  is  moderately  strong  at  the  bottom.  At  the  location,  (3.9,  3.9),  the  magnetic  field  is  close  to 
zero  at  the  bottom.  Although  the  incident  waves  cross  the  magnetic  canopy,  there  is  little  evidence 
for  mode  conversion. 

The  relative  change  in  energy  (as  a  function  of  time),  integrated  across  all  horizontal  locations, 
is  presented  in  figure  [14]  It  reinforces  the  conclusion  that  although  some  part  of  the  wave  energy 
is  transmitted  to  the  corona,  a  bulk  of  the  wave  energy  is  deposited  in  the  chromosphere  and  at 
the  base  of  the  transition  region.  This  observation  is  similar  to  the  two  dimensional  case. 

In  three  dimensions,  it  is  expected  that  Alfven  waves  are  either  excited  at  the  bottom  boundary 
or  are  generated  by  mode  conversion.  We  need  to  examine  the  numerical  results  carefully  in  order 
to  investigate  the  presence  of  Alfven  waves.  In  Carlsson  and  Bogdan  (2006),  the  authors  proposed 
that  differences  in  the  wave  structure  along  the  principal  normal  and  binormal  directions  to  the 
magnetic  field  indicates  the  presence  of  two  different  transverse  wave  modes.  This  provides  some 
indirect  evidence  for  the  existence  of  Alfven  waves.  From  figure  [l2|  we  clearly  observe  that  the 
waves  along  the  two  transverse  directions  are  different.  Does  this  indicate  the  presence  of  Alfven 
waves?  In  order  to  investigate  this  issue  further,  we  plot  the  relative  change  in  pressure  and  the 
velocity  components  parallel,  in  the  principal  normal  and  binormal  directions  to  the  magnetic 
field  in  figure  [15]  All  the  plots  are  at  time  t  =  1.4  and  we  show  the  views  from  the  top  of  the 
domain  as  well  as  from  the  side  in  figure  [l5j  It  is  well  known  that  the  pressure  should  not  change 
across  an  Alfven  wave  whereas  both  fast  and  slow  mode  waves  change  the  pressure.  A  comparison 


WELL-BALANCED  SCHEMES  FOR  WAVE  PROPAGATION  IN  NON-ISOTHERM AL  ATMOSPHERES  17 


(a)  T  =  0.9 


(b)  T  =  1.4 


(c)  T  =  1.8 


(g)  T  =  0.9 


(h)  T  =  1.4 


(i)  T  =  1.8 


Figure  12.  Results  for  the  3-D  high  frequency  planar  wave  experiment  with 
synthetic  magnetic  field  and  driving  (5.2).  Top:  Velocity  in  the  direction  of  the 
magnetic  field.  Middle:  Velocity  in  the  principal  normal  direction.  Bottom: 
Velocity  in  the  birnormal  direction  to  the  magnetic  field.  All  the  results  are  with 
a  second-order  accurate  scheme  on  a  mesh  of  400  x  400  x  800  points. 


of  the  relative  change  in  the  pressure  and  the  velocity  along  the  magnetic  field  shows  that  the 
pressure  jump  is  correlated  to  the  velocity  along  the  field  for  most  of  the  domain.  These  waves 
are  a  combination  of  fast  and  slow  modes.  In  the  region  where  there  is  no  evidence  for  waves  in 
the  direction  of  the  magnetic  field,  the  pressure  jump  seems  to  be  caused  by  a  combination  of  the 
two  transverse  wave  modes.  This  indicates  that  both  modes  are  probably  (combinations  of)  fast 
waves  and  there  is  no  direct  evidence  for  the  existence  of  coronal  Alfven  waves  in  this  particular 
case. 

5.1.2.  Transverse  driving.  We  consider  the  same  configuration  as  in  the  previous  numerical  ex¬ 
periment  but  prescribe  the  transverse  driving 

(5.3)  u  =  —  0.1(cos(27ru;t),  sin(27ro;t),  0)T,  B  =  0.1(cos(27rco’£),  sin(27ru;t),  0)T. 
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E  at  (x,y)  =  (2.0,  2.0) ,  f  =  3.0,  angle  at  tp  =  1  E  at  (x,y)  =  (2.25,  2.0) ,  f  =  3.0,  angle  at  tp  =  32 


(a)  0 ,y)  =  (2,2) 


(b)  (x,y)  =  ( 2.25,2) 


E  at  (x,y)  =  (2.2,  2.2) ,  f  =  3.0,  angle  at  tp  =  20  E  at  (x,y)  =  (3.9,  3.9) ,  f  =  3.0,  angle  at  tp  =  0 


(c)  (x,  y)  =  (2.2,  2.2) 


(d)  (x,y)  =  (3.9,  3.9) 


Figure  13.  Relative  change  in  total  energy  (as  a  function  of  height  and  time)  at 
four  different  horizontal  locations  in  the  high  frequency  planar  wave  experiment 
with  synthetic  magnetic  field  and  driving  (|5.2|).  All  the  results  are  with  a  second- 
order  accurate  scheme  on  a  mesh  of  400  x  400  x  800  points. 


E,  f  =  3.0,  angle  at  tp  =  8 


Figure  14.  Relative  change  in  total  energy,  integrated  across  all  horizontal  lo¬ 
cations,  vs.  time  in  the  3-D  vertical  driving  case. 


The  frequency  uj  =  3.0  is  considered  and  the  results  are  shown  in  figures  [T6] -  19  In  figure  [l6j 


show  waves  along  the  magnetic  field  and  in  the  principal  normal  and  binormal  directions  to  the 
magnetic  field.  As  the  forcing  is  transverse,  we  are  exciting  waves  horizontally.  In  this  case,  the 
shape  and  structure  of  the  waves  is  very  different  from  that  of  the  vertical  driving  case,  compare 
with  figure  [121  The  waves  in  the  direction  of  the  magnetic  field  have  a  spiral  motion  along  the 
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(e)  Pressure 


(f)  V(Normal)  (g)  V(Principal)  (h)  V(binormal) 


Figure  15.  Results  for  the  high  frequency  planar  wave  experiment  with  synthetic 
magnetic  field  and  driving  (5.2 ).  We  show  change  in  pressure  and  velocities  in  the 
three  orthogonal  directions  at  time  t  =  1.4.  Top:  View  down  the  top  boundary. 
Bottom:  Side  view.  All  the  results  are  with  a  second-order  accurate  scheme  on  a 
mesh  of  400  x  400  x  800  points. 


magnetic  field.  The  magnetic  field  is  at  its  strongest  at  the  center  of  the  domain  and  we  see  that 
the  waves  spiral  around  this  central  region.  On  the  other  hand,  the  waves  in  the  two  transverse 
directions  have  a  directed  upward  motion  in  the  center  of  the  domain.  As  in  the  two-dimensional 
case,  transverse  forcing  seems  to  act  like  a  localized  piston  sending  out  waves  from  the  middle  of 
the  bottom  boundary,  where  magnetic  field  is  concentrated. 

There  is  a  clear  evidence  of  mode  separation  once  the  waves  have  crossed  the  transition  layer  into 
the  corona.  As  seen  in  figure  [l6|  (at  time  t  =  1.4),  the  spiral  waves  in  the  direction  of  the  magnetic 
field  are  clearly  slow  mode  waves.  On  the  other  hand,  waves  in  the  transverse  directions  are  much 
faster.  Near  the  top  boundary,  the  fast  and  Alfven  waves  travel  at  the  same  speed.  So,  it  is  unclear 
whether  the  transverse  waves  are  fast  mode  waves  or  Alfven  waves.  We  examine  this  issue  in  figure 


17  where  the  relative  change  in  pressure  is  plotted  along  with  the  velocities  in  all  three  orthogonal 


directions.  Figure  [17|  clearly  shows  that  the  pressure  jump  is  completely  correlated  with  the  slow 
mode  wave.  On  the  other  hand,  waves  in  both  transverse  directions  do  not  have  any  pressure 
change  across  them.  This  is  a  clear  signature  of  Alfven  waves.  Furthermore,  these  transverse 
waves  are  constricted  by  the  strong  magnetic  field  which  is  another  characteristic  of  Alfven  waves 
as  fast  mode  waves  can  spread  out  across  the  magnetic  field.  Together,  these  observations  provide 
fairly  strong  direct  evidence  for  the  presence  of  Alfven  waves  in  this  configuration. 

The  energy  carried  by  waves  (as  a  function  of  time)  at  four  different  horizontal  locations  is 
shown  in  figure  [18]  There  are  marked  differences  in  the  energy  landscape  at  different  horizontal 
locations.  At  the  middle  of  the  domain  (x,y)  =  (2,2),  the  magnetic  field  is  quite  strong  at  the 
bottom  and  the  incident  waves  do  not  intersect  the  magnetic  canopy.  The  incident  wave  energy  is 
mostly  accumulated  in  the  chromosphere  and  at  the  base  of  the  transition  region.  Some  of  the  wave 
energy  is  transmitted  into  the  corona.  Significant  wave  energy  is  trapped  within  the  transition 
layer,  creating  what  appears  like  a  standing  wave.  Trailing  incident  waves  interact  with  this 
trapped  energy  and  are  unable  to  transmit  wave  energy  into  the  corona.  The  energy  landscape  at 
the  locations  (2.2,  2.2)  and  (2.25,  2)  is  quite  similar.  These  locations  have  a  background  magnetic 
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(a)  T  =  0.9 


(b)  T  =  1.4 


(c)  T  =  1.8 


(d)  T  =  0.9 


(e)  T  =  1.4 


(f)  T  =  1.8 


(g)  T  =  0.9 


(h)  T  =  1.4 


(i)  T  =  1.8 


Figure  16.  Results  for  the  three  dimensional  high  frequency  planar  wave  exper¬ 
iment  with  synthetic  magnetic  field  and  transverse  driving  (5.3).  Top:  Velocity 
in  the  direction  of  the  magnetic  field.  Middle:  Velocity  in  the  principal  normal 
direction.  Bottom:  Velocity  in  the  birnormal  direction  to  the  magnetic  field.  All 
the  results  are  with  a  second-order  accurate  scheme  on  a  mesh  of  400  x  400  x  800 
points. 


field  that  is  of  moderate  strength  at  the  bottom.  In  this  case,  there  is  evidence  of  low  energy 
waves  in  the  corona  that  are  faster  than  the  incident  waves.  We  believe  that  these  are  signatures 
of  the  Alfven  waves.  A  much  larger  proportion  of  the  wave  energy  is  transmitted  into  the  corona 
at  these  locations  compared  to  other  horizontal  locations.  At  the  location  (3.9,  3.9),  there  are  very 
small  amplitude  waves  as  the  magnetic  field  is  really  weak  at  the  bottom.  A  considerable  part  of 
the  incident  wave  energy  is  reflected  from  the  transition  layer  in  this  case.  The  relative  energy 
change  (integrated  across  all  horizontal  locations)  is  shown  in  figure  [l9j  It  shows  that  most  of  the 
energy  is  accumulated  at  the  upper  chromosphere  while  a  small  proportion  of  it  is  transmitted 
into  the  corona.  The  main  difference  with  the  vertical  driving  lies  in  the  amplitude  of  the  energy 
change  and  in  the  fact  that  a  greater  proportion  of  energy  is  transmitted  into  the  corona. 
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(e)  Pressure 


(f)  V(Normal)  (g)  V(Principal)  (h)  V(binormal) 


Figure  17.  Results  for  the  high  frequency  planar  wave  experiment  with  synthetic 
magnetic  field  and  transverse  driving  (5.3).  Top:  Velocity  in  the  direction  of  the 
magnetic  field.  Middle:  Velocity  in  the  principal  normal  direction.  Bottom: 
Velocity  in  the  birnormal  direction  to  the  magnetic  field.  All  the  results  are  with 
a  second-order  accurate  scheme  on  a  mesh  of  400  x  400  x  800  points. 


5.2.  Observed  magnetic  field.  The  above  numerical  experiments  were  performed  with  synthetic 
magnetic  fields  and  bottom  boundary  conditions.  The  real  test  of  a  code  like  SURYA  is  its 
performance  on  observed  data  sets.  For  this  purpose,  we  present  a  three-dimensional  simulation  of 
the  solar  atmosphere  with  the  same  temperature  profile  as  in  the  previous  numerical  experiments. 


The  background  magnetic  field  B  is  given  by  (2.10)  where  the  Fourier  coefficients  are  extracted 


from  the  magnetic  field  at  the  bottom  boundary.  Using  the  3d  data  model  of  Carlsson  and  Stein 
(2002),  measurements  of  the  solar  radial  magnetic  field  by  the  MDI  instrument  on  SOHO  in  1997 
are  used  to  obtain  the  magnetic  field  at  the  bottom  boundary.  This  observed  background  field  is 
depicted  in  figure  [20j  It  is  very  complex  with  a  rather  chaotic  combination  of  open  and  closed 
loop  field  lines. 


5.2.1.  Observed  driving.  We  excite  waves  at  the  bottom  boundary  by  using  a  forcing  that  was 
observed  in  the  same  part  of  the  solar  surface  by  SOHO  and  at  the  same  time  as  the  background 
magnetic  field.  The  results  of  this  simulation  with  observed  magnetic  field  and  observed  driving 
are  shown  in  figures  [21,  [22  and  23 


In  figure  [21]  we  plot  the  velocities  in  the  direction  of  the  magnetic  field  and  in  the  principal 
normal  and  binormal  directions  to  the  magnetic  field.  The  results  show  very  complicated  wave 
behavior.  There  are  strong  differences  between  waves  in  the  three  directions.  No  obvious  symme¬ 
tries  are  seen  in  the  plot.  This  could  be  attributed  to  the  nature  of  the  forcing  which  will  contain 
a  combination  of  vertical  and  transverse  driving  motions. 

It  is  simpler  to  analyze  the  energy  carried  by  waves.  We  plot  the  relative  change  in  energy  at 
four  different  horizontal  locations.  There  is  a  very  complex  spatial  variation  in  the  energy  transfer. 
At  the  location  (7.6,  2.3),  the  magnetic  field  is  very  weak.  Here  the  incident  waves  carry  energy 
up  to  the  transition  region.  Most  of  the  energy  is  dumped  at  the  base  of  the  transition  region,  and 
only  a  small  part  of  the  wave  energy  is  transmitted  into  the  corona.  At  the  location  (5.3,  8.4),  the 
magnetic  field  is  stronger  at  the  bottom  compared  to  the  previous  location.  Here,  there  is  evidence 
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E  at  (x,y)  =  (2.0,  2.0) ,  f  =  3.0,  angle  at  tp  =  1 


E  at  (x,y)  =  (2.25,  2.0) ,  f  =  3.0,  angle  at  tp  =  33 


(a)  (x,  y)  =  (2.0,  2.0) 


(b)  (x,y)  =  (2.25,2.0) 


E  at  (x,y)  =  (2.2,  2.2) ,  f  =  3.0,  angle  at  tp  =  20 


(c)  (a,  2,)  =  (2.2,  2.2) 


Figure  18.  Relative  change  in  total  energy  (as  a  function  of  height  and  time)  at 
four  different  horizontal  locations  in  the  three  dimensional  high  frequency  planar 
wave  experiment  with  synthetic  magnetic  field  and  driving  (5.3).  All  the  results 
are  with  a  second-order  accurate  scheme  on  a  mesh  of  400  x  400  x  800  points. 


Figure  19.  Relative  change  in  total  energy,  integrated  across  all  horizontal  lo¬ 
cations,  vs.  time  in  the  3-D  transverse  driving  case. 


for  fast  wave  modes  in  the  corona.  However,  these  waves  carry  a  small  amount  of  energy.  Most  of 
the  wave  energy  is  contained  in  incident  waves  that  dump  it  in  the  chromosphere.  Some  of  this 
energy  is  also  transmitted  into  the  corona.  The  magnetic  field  is  even  stronger  at  the  the  bottom 
of  the  location  (17.5,12.2).  Here,  there  is  considerable  evidence  of  fast  waves  in  the  corona.  A 
larger  proportion  of  the  wave  energy  is  transmitted  into  the  corona.  At  this  location,  there  is 
a  constant  forcing  at  the  photospheric  boundary,  but  only  a  small  part  of  the  wave  propagates 
upwards.  This  could  very  well  be  the  result  of  the  magnetic  field  that  is  almost  horizontal  at  this 
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Figure  20.  The  background  magnetic  field  for  simulations  with  observed  data. 

point.  The  behavior  at  the  location  (14.5, 16.8)  is  similar.  There  are  a  large  number  of  low  energy 
waves  in  the  corona.  There  is  energy  accumulation  at  the  base  of  the  transition  region  at  later 
times  at  this  location. 

Although  the  detailed  spatial  picture  for  energy  transfer  by  waves  is  quite  complicated,  the  gross 
(integrated  across  horizontal  locations)  change  in  energy  is  very  similar  to  that  of  the  synthetic 
magnetic  fields.  We  plot  this  integrated  relative  energy  change  (as  a  function  of  time)  in  figure [23| 
The  figure  shows  that  the  detailed  wave  dynamics  cancel  out  in  the  mean  and  the  main  features 
are  the  deposition  of  a  bulk  of  the  wave  energy  in  the  chromosphere  with  a  smaller  proportion  of 
the  energy  leaking  into  the  corona. 


6.  Conclusion 


We  simulate  waves  in  the  solar  atmosphere  using  SURYA.  This  code  is  based  on  solving  the 


equations  of  stratified  MHD  (2.1 )  using  a  recently  developed  high-order  well-balanced  finite  volume 
scheme.  The  code  uses  stable  finite  volume  schemes  resulting  from  an  upwind  discretization  of 
the  Godunov-Powell  source  terms  for  the  ideal  MHD  equations.  It  is  able  to  handle  two  and 
three  dimensional  configurations  with  diverse  background  magnetic  fields  and  driving  mechanisms, 
including  observed  data  sets. 

We  employ  SURYA  to  perform  simulations  of  waves  in  various  configurations.  In  two  space 
dimensions,  we  performed  simulations  with  a  complex  synthetic  magnetic  field  and  observe  that 

•  The  behavior  of  waves  strongly  depends  on  the  directionality  of  the  forcing.  For  vertical 
planar  forcing,  the  waves  in  the  direction  of  the  magnetic  field  are  excited  over  the  entire 
domain  whereas  transverse  waves  are  only  excited  in  regions  of  strong  magnetic  field.  On 
the  other  hand,  when  we  force  the  transverse  component,  both  normal  as  well  as  transverse 
waves  are  excited  only  in  regions  of  strong  magnetic  field.  In  the  transverse  forcing  case, 
the  driving  takes  the  form  of  a  localized  piston  at  concentrations  of  the  magnetic  field. 

•  There  is  mode  conversion  at  the  magnetic  canopy  (/3  =  1.2),  particularly  in  regions  of 
strong  magnetic  field.  Fast  waves  are  generally  excited  at  the  canopy.  Furthermore,  the 
transition  layer  serves  to  accelerate  the  waves  on  account  of  the  increase  in  sound  speed 
as  the  temperature  increases  in  the  corona.  Hence,  there  is  a  complex  co-existence  of  fast 
and  slow  mode  waves  in  the  corona. 

•  For  both  types  of  forcing,  waves  are  transmitted  into  the  corona  although  there  are  reflec¬ 
tions  from  the  base  of  the  transition  layer.  Furthermore,  waves  that  travel  horizontally 
along  the  transition  layer  are  also  excited. 
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(a)  T  =  0.9 


(b)  T  =  1.4 


(c)  T  =  1.8 


(d)  T  =  0.9 


(e)  T  =  1.4 


(f)  T  =  1.8 


(g)  T  =  0.9 


(h)  T  =  1.4 


(i)  T  =  1.8 


Figure  21.  Results  for  the  three  dimensional  experiment  with  observed  magnetic 
field  and  driving.  Top:  Velocity  in  the  direction  of  the  magnetic  field.  Middle: 
Velocity  in  the  principal  normal  direction.  Bottom:  Velocity  in  the  binormal 
direction  to  the  magnetic  field.  All  the  results  are  with  a  second-order  accurate 
scheme  on  a  mesh  of  120  x  120  x  60  points. 


There  is  a  rich  spatial  variation  in  the  energy  transferred  by  the  waves.  The  energy 
transfer  is  strongly  dependent  on  magnetic  field  strength  compared  to  the  pressure  at  a 
given  location. 

For  both  types  of  forcing,  most  of  the  wave  energy  (carried  by  incident  acoustic  waves) 
is  dumped  in  the  chromosphere  and  at  the  base  of  the  transition  layer.  Some  of  it  is 
reflected  down  the  chromosphere  and  a  smaller  proportion  is  transmitted  into  the  corona. 
The  main  difference  between  the  planar  and  transverse  forcing  seems  to  lie  in  the  fact  that 
a  greater  proportion  of  the  energy  is  transmitted  in  the  transverse  forcing  case. 

The  dynamic  behavior  of  the  waves  and  the  details  of  energy  transfer  through  waves  is 
less  dependent  on  the  forcing  frequency.  But  we  observe  that  in  the  low  frequency  case, 
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Figure  22.  Relative  change  in  total  energy  (as  a  function  of  height  and  time) 
at  four  different  horizontal  locations  in  the  three  dimensional  experiment  with 
observed  magnetic  field  and  driving.  All  the  results  are  with  a  second-order 
accurate  scheme  on  a  mesh  of  120  x  120  x  60  points. 


Figure  23.  Relative  change  in  total  energy,  integrated  across  all  horizontal  lo¬ 
cations,  vs.  time  in  the  3-D  observed  case. 


more  of  the  pressure  waves  are  stuck  in  the  transition  region  than  in  the  high  frequency 
case. 

In  three  space  dimensions,  we  perform  simulations  with  both  synthetic  and  observed  magnetic 
fields.  An  important  difference  from  the  two-dimensional  case  is  the  possibility  of  two  transverse 
wave  modes  in  three  space  dimensions.  We  follow  Carlsson  and  Bogdan  (2006)  and  plot  the 
velocity  in  the  principal  normal  and  binormal  directions  to  the  magnetic  held  in  order  to  detect 
both  transverse  wave  types.  The  three-dimensional  results  are  summarized  below. 
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•  Wave  dynamics  is  strongly  dependent  on  the  nature  of  the  forcing.  For  vertical  forcing 
at  the  bottom  boundary,  we  observe  that  waves  along  the  magnetic  field  are  present 
throughout  the  domain.  On  the  other  hand,  transverse  waves  exist  only  when  the  field  is 
strong  compared  to  the  gas  pressure.  Although  there  are  substantial  differences  between 
waves  in  the  two  transverse  directions,  we  were  unable  to  obtain  direct  evidence  for  the 
presence  of  Alfven  waves  in  the  case  of  vertical  driving.  There  are  pressure  jumps  across 
the  transverse  waves  indicating  that  they  are  fast  waves. 

•  When  the  bottom  boundary  is  forced  in  the  transverse  direction,  the  wave  dynamics  are 
very  different.  In  particular,  we  observe  a  clear  separation  between  slow  mode  waves 
and  transverse  waves  that  travel  at  a  faster  speed.  The  slow  mode  waves  spiral  along 
the  magnetic  field.  The  transverse  waves  are  directed  upward  along  the  magnetic  field. 
There  is  no  pressure  jump  across  the  transverse  waves  indicating  the  existence  of  coronal 
Alfven  waves  in  this  configuration.  This  should  be  contrasted  with  the  other  papers 
describing  numerical  simulations  like  Carlsson  and  Bogdan  (2006)  which  provide  only 
indirect  evidence  for  the  possible  presence  of  Alfven  waves. 

•  There  is  a  rich  spatial  variation  in  the  energy  transferred  by  upward  moving  waves  in  three 
space  dimensions.  The  energy  transfer  seems  to  depend  on  the  magnetic  field  strength  at 
a  given  horizontal  location. 

•  The  energy  dynamics  shows  differences  in  detail  between  the  vertical  forcing  and  the 
transverse  forcing  scenarios.  However,  the  overall  picture  is  quite  similar.  Most  of  the 
incident  energy  is  accumulated  at  the  base  of  the  transition  region  and  in  the  chromosphere. 
This  energy  is  carried  by  waves  that  move  horizontally  along  the  transition  layer.  Some  of 
the  wave  energy  is  transmitted  into  the  corona.  The  proportion  of  transmitted  to  incident 
energy  is  higher  for  the  transverse  forcing  case. 

We  also  consider  a  background  magnetic  field  and  a  bottom  forcing  that  were  observed  by  SOHO 
and  perform  a  simulation  with  this  observed  data  set  as  input.  The  wave  dynamics  are  much  more 
complicated  and  there  is  evidence  for  the  existence  of  all  wave  modes.  However,  the  gross  energy 
transfer  landscape  is  very  similar  to  the  synthetic  test  cases.  It  shows  that  a  large  proportion  of  the 
wave  energy  is  accumulated  in  the  chromosphere  providing  a  plausible  explanation  for  observed 
chromospheric  heating.  A  smaller  proportion  of  wave  energy  is  transmitted  into  the  corona.  There 
is  a  clear  evidence  for  the  existence  of  different  wave  modes  in  the  corona. 

Our  simulations  show  that  realistic  synthetic  magnetic  fields  and  synthetic  bottom  driving 
display  very  complex  wave  dynamics  and  provide  reasonable  explanations  for  phenomena  like 
chromospheric  heating,  coronal  waves,  mode  conversion  at  the  magnetic  canopy,  and  reflection  and 
refraction  at  the  transition  layer.  Furthermore,  three  dimensional  simulations  provide  evidence  for 
the  existence  of  coronal  Alfven  waves.  These  results  should  be  considered  as  a  step  towards  better 
understanding  of  solar  observations.  More  physics,  particularly  radiation,  needs  to  be  added  in 
our  code  in  order  to  provide  a  better  approximation  to  observations. 
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Appendix  A:  Numerical  schemes 


We  approximate  the  system  (2.3)  in  a  Cartesian  domain  x  =  (x,y,z)  E  [Xi,Xr\  x  [YUYr]  X 
[Zb,  Zt\  and  discretize  it  by  a  uniform  grid  in  all  directions  with  the  grid  spacing  Ax,  Ay  and  Az. 
We  set  Xi  =  Xi  -\-iAx  ,  yj  =  Yi  +jAy  and  Zk  =  Zb  +  kAz.  The  indices  are  0  <  i  <  Nx,  0  <  j  <  Ny 
and  0  <  k  <  Nz.  Set  xi+1/2  =  +  Ax/2,  yj+ 1/2  =  %  +  Ay/2  and  zk+1/2  =  zk  +  Az/2,  hence 

a  typical  cell  may  be  denoted  CiJtk  =  [a;i_1/2,  a;i+ 1/2)  x  [yj-1/2, yj+1/2)  x  [zk-1/2,  Zk+1/2) ■  The 
cell  average  of  the  unknown  state  vector  W  (approximating  U)  over  Cijik  af  time  t  is  denoted 
Wij^k{t)-  We  will  approximate  (2.1)  with  second-order  accurate  finite  volume  scheme.  Given  the 
cell  averages  W ij,k(t),  the  semi-discrete  form  of  the  second  order  scheme  is  given  by 


(6.1) 


dt 


AT  A yK 

-  1/2  -  H i^k-1/2)  +  Si,j,k  +  ^i,j,k  +  ^i,j,k  + 


The  numerical  fluxes  F,G,H  and  the  sources  S1,2,3  are  defined  below.  The  time  dependence  in 
the  above  expression  is  suppressed  for  notational  convenience. 

It  is  standard  (LeVeque  2002)  to  replace  the  piecewise  constant  approximation  W ij,k  with  a 
non-oscillatory  piecewise  linear  reconstruction  in-order  to  obtain  second-order  spatial  accuracy. 
There  are  a  variety  of  reconstructions  including  the  popular  TVD-MUSCL  limiters  (Van  Leer 
1979).  However,  we  need  (6.1)  to  preserve  a  suitable  discrete  version  of  the  steady  state  (2.7) 
and  a  standard  reconstruction  of  the  conservative  variables  does  not  lead  to  such  a  well-balanced 
scheme. 

Consequently  we  now  modify  the  novel  equilibrium  variables  based  reconstruction  algorithm  of 
Fuchs  et  al.  (2010b)  to  the  setting  of  non-isothermal  atmospheres.  Given  the  cell  averages  W ij^k 
at  any  given  time,  we  define  a  piecewise  constant  temperature  distribution  by  (2.4)  and  denote 
the  cell  temperature  as  Tij:k-  We  have  the  following  reconstruction  algorithms. 
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6.1.  Minmod  reconstruction  (MM).  Given  the  cell  values  qij,k  of  a  state  variable  q ,  denote 
the  minmod  derivative  in  the  x-direction  as 


(6.2) 


Dxq. 


minmod  ( gi+1’3''*  ,  ***- 

v  Ax 


qi  —  l,j,k 


Ax 


minmod(a,  b)  =  -(sgn(a) 


'  sgn(fe))  min(|a|,  |6|). 


The  minmod  derivatives  Dyqij ^  and  Dzqij ^  in  the  remaining  directions  are  defined  analoguously. 
Then,  a  piecewise  linear  non- os  dilatory  approximation  of  q  is  of  the  form, 

(6.3)  q  =  qi,j,k(x,  y,  z)  =  qitj,k  +  Dxq(x  -  xt)  +  Dvq(y  -  y j)  +  Dzq(z  -  zk),  (x,  y,  z)  € 

From  the  cell  values  of  p,  u  and  B,  we  define  the  minmod  slopes  by  (|6.2|)  and  obtain  the  correspond¬ 


ing  piecewise  linear  approximations  of  these  variables  by  (6.3).  However,  a  minmod  reconstruction 


of  the  pressure  does  not  lead  to  a  well  balanced  scheme.  We  need  a  novel  pressure  reconstruction, 
based  on  a  corresponding  reconstruction  of  the  temperature. 


Given  the  cell  averages  we  compute  the  cell  temperature  by  (2.4).  We  reconstruct 

the  temperature  as  a  piecewise  linear  function  taking  the  values  at  the  cell  centers.  This 


piecewise  linear  temperature  can  be  used  to  compute  the  corresponding  a  function  by  (2.6) 
However,  we  only  need  the  differences  in  a  given  by, 


(6.4)  C^,jf,/c+A  & i,j,k  — 


rZk+ A 
j  Zk 


1 


-dz  = 


Az 


■log 


(ffe+ a)\ 


a  e  {1/2,1}. 


’  zk  Ttj(z)  Ti,j,k+1  T’i,j,k  y  r,j,k  J 

Note  that  the  difference  in  a  is  always  well-defined  and  for  Tij^+i  =  Ti,j,k  degenerates  to 

(6-5)  dij^k- fA  QLi,j,k  =  (%i,j,k+ A  %i,j,k) /Ti,j,k  :  ^  £=  {1/2,  !}• 

We  use  the  a  function  to  reconstruct  the  pressure. 

As  in  Fuchs  et  al.  (2010b),  we  define 

(6.6)  L pij'k  =  log  (pijtk), 

and  compute  the  minmod  derivatives  Dx,yLp  by  (|6.2|.  A  scaled  minmod  derivative  in  the  2- 
direction  takes  the  form, 

^Pi,j,k- (-1  L Pi,j,k  L Pi,j,k  ^PiJ^k—l 


(6.7) 


Dz~Lpij,k  =  minmod 


1  &i,j,k  &i,j,k  &i,j,k  —  l 


where  the  difference  in  a  is  computed  in  (6.4).  Then  a  piecewise  linear  approximation  of  the 
pressure  is  computed  by 


(6.8) 


p(x,  y,  z)  =  pi  j  keDX1Llpi^x(kx-Xi) eDyL,P^xk(y-yj) eDzLpi,jjk(a(z)-a(zk)) , 


where  a  is  again  computed  from  (6.4).  The  cell  edge  values  of  the  conservative  variables  can  be 
easily  obtained  from  the  piecewise  linear  approximations  of  the  primitive  variables.  The  minmod 
limiter  is  one  possible  choice  among  many  reconstruction  procedures.  Other  limiters  like  the  MC, 
superbee,  ENO  and  WENO  limiters  can  be  modified  analogously 

We  denote  the  reconstructed  piecewise  linear  conservative  variables  in  the  cell  Gj,/c  as  W ij,k{%:  y , 
define  the  following  cell  edge  values, 

WiJjfe  =  ^^i,j,k{xi— 1/2?  Uj:  zk): 


W ijk  —  Wij:k(xi+ 1/2? 


Uj  :  zk  )  5 

WZ,k  =  W^(Xi,%-+i/2,^), 
wlj,k  =  w  i,jAxiiyjiZk+i/2), 


j,k(xi:  Dj— 1/2:  Zk): 
Wtj,k  =  W  ij,k{Xi:yj:Zk-l/2): 


.  ,  =  W 

vv  i,j,k  vv  * 

Tb 


and  define  the  numerical  fluxes  by 
(6.9) 

Fi+i/2,i,fe  =  F  i+i/2,j,kj  , 

Hi,j,fc+ 1/2  =  H  (w*yiiyy,w,ByW)  , 


Gi,j+i/2,k  —  G  , 
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The  value  of  the  staggered  coefficient  B  is  given  by  a  simple  evaluation, 

®(*^i+l/2  5  Vj  5  %k)  5  ®i,jf-l-l/2,fc  B(x^,  yj0k/2'>  %k)  5  -^i,j,/c+l/2  -^(*^5  Vj  5  ^+1/2)5 

This  choice  ensures  formal  second  order  accuracy  for  a  smooth  background  magnetic  field. 


6.2.  Numerical  fluxes.  Following  Fuchs  et  al.  (2010b,  2011a),  we  determine  the  numerical  flux 
Fi+x/2 ,j,k  and  the  source  term  &}  j  k  (in  the  x-direction)  from  the  (approximate)  solution  of  the 
following  Riemann  problem 


(6.10) 


Wf 


(W,BM)  =  s1(w,Bm,wA, 

W(.r.O)  =  J 

V  /  x  v  / 

l 

x  <  0, 
x  >  0, 


We  will  approximate  the  eight  waves  in  the  MHD  Riemann  problem  with  three  waves,  i.e,  two 
representing  the  outermost  fast  waves  and  a  middle  wave  approximating  the  material  contact 
discontinuity.  The  approximate  solution  and  fluxes  are  given  by 


(6.11) 


WHs  =  ^ 


WL 

if  f  <  SL, 

1 

'FL 

if  f  <  SL, 

W*L 

if  sl  <  f  <  sm, 

03 

II 

n 

if  sl  <  f  <  sm, 

W  *R 

if  sM  <  f  <  sR, 

n 

if  SM  <  f  <  SR 

if  sR  <  f , 

[Fr 

if  Sr  <  J. 

We  set  7Ti  =  p-\-  B‘2^2Bs .  The  outer  wave  speeds  sl  and  sR  model  the  fast  magneto-sonic  waves  and 
are  defined  as  in  Gurski  (2004).  In  order  to  describe  the  solver,  we  need  to  determine  the  speed 
of  the  middle  wave  %  and  the  intermediate  states  W£,  W^.  The  middle  wave  models  a  material 
contact  discontinuity.  Hence,  the  velocity  field  and  the  tangential  magnetic  fields  are  assumed  to 
be  constant  across  the  middle  wave.  This  allows  us  to  define  u*  =  u^  =  uj^,  H2  =  B2L  =  B2R 
and  H3  =  B3L  =  B3R.  The  normal  magnetic  field  Hi  is  not  assumed  to  be  constant  but  jumps 
only  across  the  middle  wave  (modeling  the  linear  degenerate  “divergence  wave”  and  Hi  is  constant 
across  the  outer  waves.  The  intermediate  states  are  determined  by  local  conservation  across  the 
two  outermost  waves  and  the  middle  wave  resulting  in, 


(6.12) 


s<Tw:  -  F;  =  -  F„,  smW*r  -  sMW*L  =  F*r-F 


.1,* 


where  a  E  {H,  R}  and 


(6.13) 


/  0  \ 

bIr-b\l\ 

2  t 

~  {£>2 j  (Bir  ~  Bil) 

-  (^3)  {Bir  ~  Bil) 

_u*  (^lR  ~  Bil ) 

—  -  (u*2B*2  +  u*3B*3)  (B1R-B1L)J 


This  amounts  to  integrating  the  source  s1  across  the  wave  fan. 
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Applying  the  conservation  relations,  we  obtain  (check  Fuchs  et  al.  (2010b),  section  3.1.2  for 
details)  the  following  intermediate  states, 


Po  =  Pe— - 7r*0  =  tti e  +  Po(u le  ~  sq)(uiq  -  sM ),  0  G  {A,  i?}, 


sm  =  u  i  = 


%  —  ^0 

TTlR  -  TTlL  +  />R^1R(^1R  “  5fl)  ~  /AL^1l(^1L  “  «l) 


Pi?(^lR  -  Sr)  ~  Pl{uiL  -  SL) 

a  G  {2,3} 


*  C^cr  —  *  ada  £c( 

—  r  ITT  i  -l)^.  — 


—2 


«C  +  C2  ’  «C  +  C2 

1  /  D~ 

El  =  - (  ^<9  (^10  50)  +  TTiqUio  -  TtIqSm  H - ^  (^10  -  5m) 

5m  —  50  V  2 

+  {E\q)  (^B 20^29  +  B30U30  -  B2eu*2e  -  B3eu%e )) ,  0  G  {L,  i?}, 


ca  =  PrUcjR  (uir  —  Sr)  —  pLUaL  {uil  —  Sl)  —  (BiRBaR  —  BiLBaL)  , 

d(j  =  (i^ii?  —  Si?)  —  BaL  (uil  —  5l)  —  ( BilUcjL  —  BiRUaR)  , 
a  =  PR  (uir  sR )  pL  (u1L  -  Sl)  ,  C  =  5r  -  5l,  £  =  -  B1L. 


The  intermediate  fluxes  are  obtained  in  terms  of  the  intermediate  states  by  local  conservation 

(p2), 

El  =  Fl  +  sL(Wl  -  WL),  =  FR  +  sr(W*r  -  W R). 

The  discrete  source  term  takes  the  form, 

(6’14)  Si’"  =  Si-l/2X(SM,i-l/2>0)  +  Skl/2^(SM,i  +  l/2<0)5 


where  ,2  is  defined  in  analogy  to  (|6.13 ). 

For  the  three  dimensional  form  of  the  equations,  the  fluxes  G,H  and  the  sources  S2  and  S3 
can  be  defined  analogously.  The  discretized  gravity  source  term  S9  is  given  by, 


(6.15) 


S9  ■ , 

l,],k 


0,0, 


pUk 


■pbi,j,k 


A.z 


,0,0, 0,0  ,-fijik(v%)u,kg\ 


6.2.1.  Boundary  conditions  for  the  second  order  scheme.  The  boundary  is  treated  in  the  following 
way.  We  need  to  specify  two  layers  of  ghost  cells  in  each  direction  for  a  second  order  scheme.  We 
use  periodic  boundary  conditions  in  the  x-  and  y-  directions,  i.e.,  for  1  <  j  <  Ny  and  1  <  k  <  Nz, 
we  have 
(6.16) 

Wo,j,fc  =  W  Nx,j,k]  W-l,j,k  =  Wjva.+l,j,fe  =  Wl  J,fe,  WNx+2,j,k  =  W2  }j,jfe 

The  ghost  cell  values  in  the  ^/-direction  can  be  defined  analogously. 

In  the  ^-direction,  we  use  second-order  extrapolated  Neumann  boundary  conditions  for  the 
velocity  and  the  magnetic  field,  i.e.,  for  w  =  {u,  B}, 


(6.17) 


w  i,j,d  —  ™i,j,Nz+2+d  —  wi,j,Nz 


for  1  <  i  <  Nx,  1  <  j  <  Ny  and  d  G  {0,-1}  in  order  to  define  the  values  in  the  ghost  cells. 

The  pressure  and  the  density  in  the  ghost  cells  are  extrapolated  in  terms  of  its  logarithm 
L p  =  log(p)  and  L p  =  log(p)  according  to  (|6.8|)  and  simplify  to 


(6.18) 


Pi,j,d  =  Pi,j,ie{ai’i’d  ai^H, 
Pi,j,d  =  Pi,j,ie^’d~ai’^l)/H, 


Pi,j,Nz+2+d  =  Pi,j,Nze.{°‘i’i’N*+2+d  ai’ 
Pi,j,Nz+2+d  =  Pi,j,Nze('°‘i’i-N*+2+d-ai’ 


,1,1 vJ/H 

1 

,1,nz)/H 

5 


where  the  differences  in  a  are  given  by  (6.5).  This  amounts  to  using  a  scaled  version  of  the 
extrapolated  Neumann  type  boundary  conditions  of  Fuchs  et  al.  (2010a)  for  the  primitive  variables. 
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6.2.2.  Time  Stepping.  The  standard  scheme  for  a  first  order  approximation  in  time  is  the  forward 
Euler  time  stepping,  formally  written  as 


W 


n+ 1 


=  W] 


,3,k 


A  tnT 


n 

i,j,k 


where  is  defined  in  (|6.1|).  For  second-order  schemes,  we  use  the  second-order  strong-stability 

preserving  Runge-Kutta  (SSP)  time  stepping  (Gottlieb  et  al.  2001) 


W 

W 

w 


*  =  Wn  -L  A  tn  'Fn 

i,j,k  i,j,k  '  ^  i,j,k"> 


**  .  =  w*  •  u 

h3,k  i,3,k 


+  A  tnT* 


i,j,ki 


n+i  =  -fWn  •  , 

i,j,k  2V 


The  time  step  is  determined  by  a  standard  CFL  condition. 

The  scheme  (6.1 )  as  constructed  above  is  shown  to  be  second-order  accurate  and  is  well  balanced 
i.e,  it  preserves  a  discrete  version  of  the  steady  state  (2.7).  The  proof  of  these  properties  is 
presented  in  Fuchs  et  al.  (2011b). 
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